Interacting electrons on trilayer honeycomb lattices 
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Few-layer graphene systems come in various stacking orders. Considering tight-binding models 
for electrons on stacked honeycomb layers, this gives rise to a variety of low-energy band structures 
near the charge neutrality point. Depending on the stacking order these band structures enhance 
or reduce the role of electron-electron interactions. Here, we investigate the instabilities of inter- 
acting electrons on honeycomb multilayers with a focus on trilayers with ABA and ABC stackings 
theoretically by means of the functional renormalization group. We find different types of compet- 
ing instabilities and identify the leading ordering tendencies in the different regions of the phase 
diagram for a range of local and non-local short-ranged interactions. The dominant instabilities 
turn out to be toward an antiferromagnetic spin-density wave (SDW), a charge density wave and 
toward quantum spin Hall (QSH) order. Ab-initio values for the interaction parameters put the 
systems at the border between SDW and QSH regimes. Furthermore, we discuss the energy scales 
for the interaction-induced gaps of this model study and put them into context with the scales for 
single-layer and Bernal-stacked bilayer honeycomb lattices. This yields a comprehensive picture of 
the possible interaction-induced ground states of few-layer graphene. 

PACS numbers: 71.27. +a, 71. 10.Fd,71.30.+h,73.21.Ac,75.70.Cn 



I. INTRODUCTION 

In the last years, the electronic properties of few- 
layer graphene systems have been studied extensively-^. 
Experiments on graphene systems with two 2 - - — and 
three layers^— aimed to clarify the role of many-body 
interaction ;) 18 ' 19 . These are expected to be more effective 
due to the enhanced low-energy density of states and thus 
correlated electronic ground states may be realized. For 
bilayer graphene the experimental studies agree on the 
correlated nature of the ground state, but disagree on 
the symmetries of the underlying order, a topic that is 
matter of current debate. At temperatures below 5K, 
transport experiments found gap openings of about 2- 
3 meV— — . For trilayer graphene the stacking order is cru- 
cial for the electronic properties^— . In Ref. [jj], it was 
found that neutral trilayer graphene with ABC stacking 
shows many-body correlations with a pronounced gap of 
~ 6 meV while graphene trilayers with ABA stacking do 
not show a gap. In view of this experimental situation a 
better theoretical understanding of the many-body insta- 
bilities and the nature of the correlated phases which are 
candidates for possible trilayer graphene ground states is 
required. 

Here, we employ a functional renormalization group 
(£RG) approach (for a recent review, see Ref. H3) to ad- 
dress the problem of competing instabilities on honey- 
comb trilayers in an unbiased way. We explore a region 
of the phase diagram with a range of interaction param- 
eters with density-density repulsions up to the second 
nearest neighbor that are motivated by the ab-initio val- 
ues as proposed in Ref. |3l|. In a previous study, our 
RG approach was applied to the AB stacked honeycomb 
bilayer— to investigate the phase diagram as a function 



of local and non-local interaction parameters and com- 
plemented by a more quantitative study for pure onsite 
interactions combining mean-field theory, functional RG 
and quantum Monte Carlo (QMC) techniques^ 3 -. These 
studies clarified the relevance of various phases on the bi- 
layer that have been subject to previous theoretical stud- 
ies by different authors 3 ^—, e.g. the antiferromagnetic 
spin density wave (AF-SDW), two kinds of charge den- 
sity waves (CDW, CDW 3 ) and a quantum spin Hall state 
(QSH), depending on the model parameters when the full 
four-band model with the available ab-initio estimates 
for the model parameters is used. The combination with 
QMC simulations^ 3 - showed the compatibility of fRG and 
QMC methods for the investigated range of parameters, 
and allowed to gauge the quantitative precision of the 
fRG data at least partially. 

The instability toward an interaction-induced quantum 
spin Hall (QSH) state is of particular interest in con- 
nection with realizing topological insulators in graphene. 
While the original proposal by Kane and Mele^, extend- 
ing earlier work by Haldano^ with a mass term due to 
spin-orbit interactions turned out to too small gaps, it 
was later argued^ 8 - that second-nearest neighbor repul- 
sions could lead to a interaction-induced mean-field hav- 
ing the same effect as the Kane-Mele mass term. In 
the single graphene layer, however, a nonzero interac- 
tion strength is needed to open any kind of gap. So 
far there is no experimental evidence for this. In the 
bilayer system, the same instability occurs for arbitrar- 
ily small interactions of appropriate distance dependence 
(see e.g. Ref. HI), and one basically gets two copies of the 
Kane-Mele QSH state coupled by the interlayer terms. 
The sign of the order parameter gives rise to distinct 
choices with interestingly different properties^. While 
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this is an interesting many-body state, the correspond- 
ing spin-polarized edge modes would not be topologically 
protected in a strict sense as the bilayer-doubling now 
permits time-reversal invariant single-particle terms that 
would localize the edge states^. However, if one finds the 
same state in the trilayer system, one would again have 
an odd number of Kramers pairs at the edges or edge 
states per spin, and the topological protection would keep 
alive at least one edge mode. In this work we show that 
the fRG suggests that the trilayer QSH state is not un- 
likely for realistic parameters, and that it can even occur 
at sizable energy scales ~ 10 meV. 

This paper is organized as follows. In Sec. [Ill we in- 
troduce simple tight-binding Hamiltonians for ABC and 
ABA trilayer honeycomb lattices and discuss the result- 
ing band structures depending on the stacking order. 
Further, we introduce the interacting part of the Hamil- 
tonian with onsite, nearest and next-to-nearest neigh- 
bor interaction terms. In Sec. IIIII we describe the fRG 
method and discuss the patching scheme as well as the 
employed approximations, together with the relations to 
other methods. In Scc. lIVI we analyze the emergent effec- 
tive interactions from the fRG flow. We classify the in- 
stabilities and ordering tendencies in ABA and ABC tri- 
layer honeycomb lattices for interacting electrons. With 
this information, we can construct the tentative phase 
diagram from fRG in the space of interaction parame- 
ters. We finish with a discussion of the energy scales for 
the gaps in the electronic spectrum due to the ordering 
phenomena in Sec. HVDl and some conclusions in Sec. fVl 



II. MODEL HAMILTONIANS 




FIG. 1: (a) Sketch of ABC stacked trilayer graphene with 
trigonal warping terms, (b) Sketch of ABA stacked trilayer 
graphene with trigonal warping terms. The hoppings are 
shown in side-view. 



The index s =t, 4- specifics the electron spin. With these 
preliminaries we can write down the single-layer tight- 
binding Hamiltonians in the first two layers, 



~7o £ 



C b 1 ,s,R Ca i, s ,R+f>i 



C b 2 ,s,R-8i Ca 2,s,R 



h.C. 



(1) 

(2) 



s.R.Si 



and a stacking dependent Hamiltonian for the third layer, 
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where we have attributed an additional index c to H 3c , 
with c = for the ABA stacking and c = 1 for the ABC 
or chiral stacking. In graphene and few-layer graphene 
systems the hopping 70 = t « 3 cV, sec e.g. Refs. I40ll5ll . 
Next we introduce the intcrlaycr hoppings, 71, between 
sites that lie on top of each other and connect adjacent 
layers, 



A. Basic Hamiltonians 

Here, we construct the non-interacting part of the 
tight-binding Hamiltonians for ABC and ABA stacked 
trilayer honeycomb lattices. The position vectors of the 
two-dimensional bipartite lattice structure are called R, 
with each R having three nearest neighbors at positions 
that can be characterized by the three vectors 8 n with n € 
{1, 2, 3}. Explicitly, the 5 n are given by Si = ^§^e x +^e y , 



Here, a is the 



S 2 = — + §e y and <5 3 - 

distance between two neighboring lattice sites and the 
vectors point from the B-sublattice to the A-sublattice. 
The tight-binding Hamiltonian for ABC/ ABA stacked 
trilayers is composed out of single layer Hamiltonians for 
the in-plane hoppings, perpendicular hoppings between 
different layers and remote hoppings between the layers 
including a certain planar distance. We denote the differ- 
ent lattice sites by Ai, Bi, A 2 , B 2 , A 3 , B 3 , where A and 
B specify the sublattice and the index numbers the layer 
as shown in Fig. [TJ Accordingly, we define the annihi- 
lation operators c G s f and the creation operators - 
of an electron at position f and the layer and sublattice 
are encoded in the subscript o 6 {ai, bx, a 2 , b 2 , 13, ^3}- 
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where again we accounted for the different stackings by 
introducing H 23 c . Ab initio values for 71 = t± are avail- 
able for graphite2&£2, t± w 0.4 eV , and ABC trilayer 
graphene^, t± w 0.5 eV. In Fig. [TJ we also show the 
more remote hoppings 72,73,74,75 whose effect will be 
discussed in the following section. In the present model 
study we will ignore these terms for most explicit cal- 
culations. For the discussion of the energy bands, we 
introduce the Fourier transform 
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where Af is the number of unit cells and k is an element 
of the first Brillouin zone (BZ). We can write the tight- 
binding Hamiltonians in this approximation in Fourier 
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FIG. 2: Sketch of ABC and ABA stacked trilayer graphene 
band structure with t = 3eV, t± = 0.4 eV. Top left panel: 
ABC trilayer dispersion with k y = 0. Top right panel: ABC 
trilayer dispersion close to K, K' with k y — 0. Bottom left 
panel: ABA bands close to K, K' with k y = 0. Bottom right 
panel: ABA bands close to K, K' with k y = with remote 
hoppings 72,73,74,75 and onsite energy 5 and the numerical 
values as given in the main text. 



space as 
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where dk = X)n gx p(^ • S n ). The resulting energy bands 
are depicted in Fig. [2] The dispersion for ABC trilayers 
is very flat near the Fermi level with cubic wavevector de- 
pendence close to the K and K' points. Hence one should 
expect an enhanced role of interactions as compared to 
ABA trilayers and AB bilayers where the wavevector de- 
pendence is quadratic. 



B. Remote hoppings in trilayers 

In this section, we want to discuss the effect of remote 
hopping terms on the band dispersion and possible im- 
plications for instabilities. For the ABC trilayer system 
the density of states close to the Fermi energy has a van 
Hove singularity ~ e -1 / 3 due to the cubic band cross- 
ing point. This enhances the role of interaction effects, 



leading to high critical scales for ordering tendencies. In 
our model study, we therefore choose to ignore remote 
hoppings that change the topology of the bands below 
an energy of order 10 meV— , which is of the order of the 
measured gap in the experiment. Of course, it is an in- 
teresting question whether these terms affect the nature 
of the instability in trilayer graphene. However, we find 
in our computations that instabilities occur already on 
an energy scale that is higher than or at least compara- 
ble to 0.004£ ~ 10 meV. This serves as an a posteriori 
justification for dropping the remote hopping terms and 
proves the consistency of our approach. 

In ABA trilayer the situation is different. Here, the 
band structure close to the K-, A''-points is separated 
into a J = 1 (linear) and a J = 2 (quadratic) sub- 
band. The remote hoppings induce separate gaps for 
these subbands, however, with an individual energy shift 
of the subbands. This destroys the particle-hole nest- 
ing and the associated instabilities are suppressed (at 
least for small interaction terms). Explicitly, the tight- 
binding Hamiltonian in presence of the important re- 
mote hoppings2&22£2, 72 = -0.02 cV, 73 = 0.3 eV, 74 = 
0.04 eV, 75 = 0.04 eV, 5 = 0.05 eV, can be written as 



H* 



(10) 

Within the fRG approach, we will show that also for this 
system flows toward many-body instabilities occur, how- 
ever, only beyond a critical interaction strength which, 
similarly to the single-layer case, most probably is larger 
than the interaction strength in the real material. 



C. Interaction terms 

In order to investigate the instabilities that are possible 
on the trilayer honeycomb lattice for interacting electrons 
we will take into account a number of different interaction 
terms, most importantly an onsite or Hubbard interac- 
tion U, a nearest neighbor intralayer interaction V\ and 
a next-to-nearest neighbor intralayer interaction Vi. For 
these interaction parameters ab initio parameters from 
constrained random phase approximation (cRPA) com- 
putations are available^ and we take those values as a 
motivation for the investigated range in the phase dia- 
gram. The interaction Hamiltonian reads 

Hi = U nj^n^x + Vi }^ n itS nj tS i (11) 

+ V 2 ^ n i,s n j,s> , 
((i,j)),s,s' 

where n^ s = c\ s jC , s ,i and i,j run over the lattice sites, 
but pairs are included only once. The unitary transfer- 
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mation from the orbital to the band degrees of freedom 
diagonalizing Hq L or Hj h is performed numerically and 
has to be carried out on Hi as well. This adds 'orbital 
makeup' to the interaction terms in band representation, 
leading to an additional angular dependence of the inter- 
actions near the K, AT'-points. In order to resolve this 
dependence we use an angular patching of the interaction 
terms as explained in Sec. IIIII 



III. FRG METHOD 

We employ a functional rcnormalization group (fRG) 
approach for the one-particle-irrcducible vertices with 
a momentum cutoff (for a recent review on the fRG 
method, see^). In this scheme, an infrared regulator 
with energy scale A is introduced into the bare propaga- 
tor function, 



G (oj,k,b)^G A (u;,k,b) 



C A [e{k,b)\ 
icu — e(k, b) 



(12) 



Here, ui is the Matsubara frequency, k the wavevector, 
b the band index and e(k, b) is the single-particle disper- 
sion. As the spin is conserved the free propagator is spin 
independent. The cutoff function is chosen to enforce an 
energy cutoff, which regularizes the free Green's function 
by suppressing the modes with band energy below the 
scale A, 



C A [e(M)]«e (\e(k,b)\-A 



(13) 



For better numerical feasibility the step function is 
slightly softened in the actual implementation. With 
this modified scale-dependent propagator, we can define 
the scale-dependent effective action Ta as the Legendre 
transform of the generating functional Q\ for correlation 
functions, cf. Ref. [5^. The RG flow of T\ is generated 
upon variation of A. By integrating the flow down from 
an initial scale Ao, which is in our case chosen as the 
maximum energy of all bands, to the infrared A — > 0, 
one smoothly interpolates between the bare action of 
the system and the effective action at low energy. In 
order to limit the numerical effort we employ the fol- 
lowing approximations. First, the hierarchy of flowing 
vertex functions is truncated after the four-point (two- 
particle interaction) vertex. Secondly we neglect the fre- 
quency dependence, by setting all external frequencies 
to zero, as we are interested in ground-state properties. 
The general coupling function, which depends on three 
momentum and Matsubara frequency indices (the fourth 
index is fixed by conservation) is thus replaced by the 
coupling function V A (fa, fa, fa, 64) in band representa- 
tion with ki = (ki,bi) or in orbital representation by 

V A (k\, fa, fa, 04) with ki = (ki,Oi). Third, we neglect 
self-energy corrections. This approximate fRG scheme 
then amounts to an infinite-order summation of one-loop 
particle-particle and particle-hole terms of second order 



in the effective interactions. It allows for an unbiased 
investigation of the competition between various correla- 
tions, by analyzing the components of V A (ki, fa, fa, 04) 
that create instabilities by growing large at a critical scale 
hj&. With the approximations mentioned above, this 
procedure is well-controlled for small interactions. At 
intermediate interaction strengths we still expect to ob- 
tain reasonable results. In any case, the fRG takes into 
account effects beyond mean-field and random phase ap- 
proximations. 

The fRG calculation is performed in the band basis 
in which the free part of the Hamiltonian is diagonal. 
The free Hamiltonian can be diagonalized by a unitary 
transformation of the form 
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where the index o includes all six sublattices in the 
three layers and the index b denotes the correspond- 
ing bands. In the interaction part of the Hamiltonian 
the coupling function V A {k\, fa, fa, 04) in orbital space 
is chosen so that it reproduces the interaction Hamilto- 
nian of Eq. (fTTj) . For application of the fRG in the band 
basis we have to transform the interaction via 
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Thus, in the band basis the interaction ver- 
tex V{fa, fa, ^3, ^4) acquires a pronounced mo- 
mentum dependence due to the extra prcfactor 
u t r u , u* - u* - , which is sometimes also 

o 1 b 1 ,k 1 o 2 b 2 ,k 2 o 3 b 3 ,k 3 o 4 b4,fc 4 ' 

referred to as 'orbital makeup'. Previous studies have 
shown that this has a crucial impact on the fRG flow by 
allowing for unconventional instabilities^. 

With the approximations discussed above the fRG flow 
for the coupling function reads 

±V A {fa,fa, fa,b 4 ) = T A p + T A H4 + T A HiCr , (17) 

with the particle-particle channel 
Tpp{k\,fa,fa,b A ) = -—— \ dkS^\ 

VbzJ Tt' 

V A {k 1 ,fa,k,b r )L A {k,qpp)V A (k,q P p,faM) \ , (18) 
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FIG. 3: Left panel: Interaction vertex labeled with the spin 
convention (upper diagram). Below, the loop contributions 
to the flow of the interaction vertex including the particle- 
particle-diagram (a), the crossed particle-hole-diagram (b) 
and the direct particle-hole diagrams (c). Right panel: Patch- 
ing of the Brillouin zone as explained in the text. 

the direct particle-hole channel 

Tp H .d( k i, k 2,k 3 ,h) = -— ^ J dk^2[ 

b,b' 

- 2V A (k u k, k 3 , b')L A (k, qPH,d)V A (qpHA, *2, K 64) 
+ V A (k, fci, fc 3 , b')L A (k, qp H ,d)V A {qPH4, k2,k, 64) 
+ V A (k ll k, fc 3 , b')L A (k, q PH . d )V A (k 2 , q PH . dl k, 64)] , 

(19) 

and the crossed particle-hole channel 
Tp H ^ cr {kuk 2 ,k 3 ,b 4 ) = -— |- J dk J2[ 

V A {k, k 2 ,k 3 , n')L A (k, qp H , C r)V A {k u qpH,cr, k, 64)] , 

(20) 

where k = (k,b) collects the wave vector and the band 
index. As mentioned above, external frequencies are set 
to zero loi = lu 2 = lo 3 = uj 4 = 0. qpp = —k + k\ + 
£2, qpH.d = k + fci - fc 3 , qpH.cr = k + k 2 - k 3 arc the 
wavevectors of the second loop line. The band index of 
the second loop line is denoted with b' . The frequency of 
the second line is fixed by frequency conservation to be 
— uj in the particle-particle diagram and uj in the particle- 
hole diagrams. Vbz 1st the volume of the BZ. The loop 
kernel is given by 

£ A (M') = ^ [G£(fc)G£(fc')] , (21) 

where in our approximation self-energy corrections are 
neglected, i.e., the full propagator is identical to the free 
propagator. 

The wavevector dependence of the interaction vertex 
is simplified by discretization. The BZ is divided into N 
patches with constant wavevector dependence within one 



patch, so that the coupling function has to be calculated 
for only one representative momentum in each patch. 
The representative momenta for the patches are chosen 
to lie close to the Fermi level. The patching scheme is 
shown in Fig. [3l with N = 24. Each of the four mo- 
menta in V A (ki, k 2 ,k 3 , 64) is additionally equipped with 
a band index. Momentum conservation fixes one of the 
four wavevectors. Altogether this results in a 6 4 x N 3 
component coupling function V A . 

We start the fRG flow at the initial scale Ao which is in 
our case chosen as the maximum energy of all bands. We 
then integrate out all modes of these bands by decreas- 
ing A. In typical flows some components of the effective 
interaction vertex become large and diverge at a critical 
scale A c > 0. In this work we use the scale at which the 
interaction vertex exceeds a value of the order of 10 times 
the bandwidth as an estimate for the critical scale. The 
precise choice of this value has only a minor effect on the 
extracted critical scale, as the couplings grow very fast 
in the vicinity of the divergence. 

The divergence is strictly speaking a (physically mean- 
ingful) artifact caused by the neglect of the self-energy in 
the flow. With self-energy correction a gap would open 
up or some other modification of the low-energy spectrum 
would take place, and the flow would be regularized. This 
is all well known from the Cooper instability in super- 
conductors. Our analysis here tells us in which channel 
ordering occurs most prominently. The pronounced mo- 
mentum structure of the vertex near the critical scale 
can be used to extract an effective Hamiltonian for the 
low-energy degrees of freedom. This is used to determine 
the leading order parameter of a given instability. Fur- 
thermore, the scale A c can be interpreted as an estimate 
for ordering temperatures, if ordering is allowed by the 
Mermin- Wagner theorem, or at least as the temperature 
below which the dominant correlations should be clearly 
observable. Furthermore, one can understand A c as en- 
ergy scale for the modification of the spectrum, typically 
by a gap. 

In this work we study the flow at temperature T = 0. 
We find flows to strong coupling with non-zero criti- 
cal scales A c for all choices of non- vanishing interaction 
terms provided there is a non- vanishing density of states 
at the Fermi level of the coupled layers. 



IV. INSTABILITIES AND PHASE DIAGRAM 

A. ABA and ABC trilayer Hubbard model 

Let us start the description of the fRG results with the 
case of onsite interactions only, i.e. U > 0, V\ = V2 = 0. 
We limit the study to the charge-neutrality point, i.e. 
with Fermi points at K and K' in the Brillouin zone. 
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FIG. 4: Effective interaction vertex near the critical scale in 
the AF regime in units of t. Left Panel: Orbital combina- 
tions with Oi = 02 = 03 = 04. The numbers on the axis 
specify the number of the patch as shown in Fig. [3] On the 
horizontal axis the wavevector ki can be read off and on the 
vertical axis we enumerate fea • &3 is fixed on the first patch, £4 
then follows from momentum conservation. Here, we see that 
the sharp vertical structure (fci = £3) comes with the dou- 
ble magnitude as the horizontal structure (&2 = kz). Middle 
Panel: Effective vertex function for the orbital combination, 
where 01 = 03,02 =04/01 and if 01 G {01,02,03} then 
02 G {61,62,63}. Right panel: Effective vertex function for 
the orbital combination, where 01 = 04,02 = 03 7^ Oi and if 
01 G {01,02,03} then 02 £ {61,62,63}. If in the second and 
third plot both 01,02 £ {01,02,03} or € {61,62,63} the sign 
of the vertices changes. 



1. Simplified band structures 

First, we want to investigate the simpler band struc- 
tures when all remote hopping terms are neglected and 
only t ^= and t± ^ 0. Then, running the fRG flow 
with pure onsite interaction U 7^ in the ABC and 
ABA stacked trilayer, we observe an instability toward 
an antiferromagnetic spin density wave (AF-SDW) with 
a typical signature of the interaction vertex near the in- 
stability as shown in Fig. @] The leading part of effective 
interaction corresponding to the clearly discernible sharp 
structures in wave vector space reads in this case 



0,0' 



t o to'^o = °q=0 



(22) 



-c , ,-. The 

o,s,k o,s',k 



withVoo* >0 a ndS| =0 = iEg, s y 
fact that the above Hamiltonian only contains the q = 
component means that the effective interaction has be- 
come infinitely-ranged^. The factors e G depend on the 
orbital, e G = +1 for o G {01,02,03} and e = —1 for 
o G {o 1; o 2 , 03}- This sign structure implements the stag- 
gering of the interaction within the unit cell appropri- 
ate for antiferromagnetic interactions. Note that this 
paramctrization holds in both cases, for the ABC as well 
as for the ABA stacking. 

A mean-field decoupling of Haf results in an AF spin 
alignment in each layer where a net spin (e.g. 'up') mo- 
ment is located on the Ai-, A2- and A3 sublattices, and 
an opposite net spin ('down') moment on the Bi-, B2- 
and B 3 -sublattices. 

The critical scale A c function of the onsite inter- 
action U for ABC and ABA honeycomb trilayers with 
model hopping parameters t = t± is shown in Fig. [5] to- 
gether with the critical scales of single- and bilayer hon- 
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FIG. 5: fRG critical scale for the singlelayer (dotted purple), 
the AB bilayer (dot dashed blue), the ABA trilayer (dashed 
orange) and the ABC trilayer (solid red). These results have 
been obtained with t± — t and all higher hopping terms set 
to zero. The inset shows the same data in a semilog plot over 
t/U. 



eycomb lattices with the same hopping parameters. This 
choice of band parameters takes us beyond the regime of 
realistic parameters for TLG, but pronounces the charac- 
teristic features of the ABC and the ABA stacking close 
to the K, K' points and therefore allows to study the dif- 
ferences of the various honeycomb stacks more explicitly. 
Furthermore, it allows to compare to recent QMC results 
for the bilayer system^. We add a systematic study of 
the dependence on t± below. Most importantly, we ob- 
serve that in the ABC trilayer, the critical scale decreases 
more slowly as compared to AB bilayer when the onsite 
interaction is decreased. While in the case of ABA tri- 
layer and AB bilayer with quadratic band crossing points 
the functional dependence of A C (U) can be fitted by an 
exponential decay ~ exp(— a/U) (cf. inset of Fig. [5]), 
this does not hold for the ABC trilayer case. Instead, 



-1/3 



one 



at small U, based on the density of states ~ e 
would naively expect a behavior A c ~ U 3 . This is how- 
ever not reproduced by our data, presumably due to the 
influence of the high energy sector, i.e. additional bands. 
We expect that the leading A c ~ U 3 dependence might 
be recovered at smaller U and thus smaller A, which is 
numerically hard to access due to the rapidly decreas- 
ing critical scale. From this analysis we conclude that 
in ABC trilayers an interaction-induced gap, which we 
expect to be of the order of the critical scale A c , may 
be considerably larger than in the other structures. This 
also implies a more stable correlated ground state. All 
this seems compatible with recent experiments^. 



2. Inclusion of remote hoppings in ABA trilayers 

When the remote hoppings 72,73,74,75 and the on- 
site energy 5 in ABA trilayers are taken into ac- 
count, the band dispersion is deformed considerably 
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and the particle-hole nesting of the band structure is 
destroye d 26 ' 27 , see Fig. [5] In this case, one should ex- 
pect that in the weakly interacting limit the tendency 
toward instabilities will be strongly reduced. However, 
it is interesting to study what happens in the case of 
intermediate and larger onsite interactions. 

An accurate treatment of the dispersion with remote 
hoppings at lowest scales would require the implemen- 
tation of a new patching scheme to resolve the vicinity 
of the K and K' points, i.e. the fact that the Fermi 
surface is now not restricted to a single point. This is 
beyond the scope of this work. Therefore we now stop 
the flow at the energy scale A* at which the bands be- 
come non-monotonous. This procedure is routinely done 
in parquet- and g-ology studies of imperfectly nested 
bands^. Of course it leaves the scales below A* unin- 
tegrated, but as the dispersion at these lowest scales is 
not nested and partially gapped, we do not expect signif- 
icant contributions of these modes to the flow. Hence we 
expect that the so-obtained estimate for critical scales 
and U c is already qui te goo d. For the hopping pa- 
rameters given in Refs. [2611271 the energy scale at which 
the dispersion is not monotonous any more is given by 
A* - 10meV~ 0.004*. 

However, it is important to notice that the remote hop- 
pings already change the dispersion far above the scale 
A*. This explains why for the ABA stacking (here from 
now on called ABA* stacking when the remote hopping 
are included) , the effect of the remote hoppings is quite 
drastic. Studying again the case of onsite interactions U 
only, we observe clearly diverging susceptibilities only for 
onsite interactions U > 2.6t at critical scales well above 
0.004i. For smaller interactions, the couplings grow very 
slowly and no divergences at finite scales above A* can 
be identified. In Fig. [5] we show the fRG results for the 
critical scales for the ABA* trilayer with remote hopping 
terms. For comparison, we also show the curves for the 
single- and bilaycr system. As the similarity to the sin- 
gle layer is strong, this analysis suggests a critical onsite 
interaction U c ~ 2.6i above which a many-body instabil- 
ity can occur in the ABA* stacking with remote hopping 
included. We would like to add that the fRG in the 
present approximation has the tendency to overestimate 
critical scales. Therefore we would expect the true U c to 
be slightly larger. For instance, in the single layer system 
for onsite interactions only, QMC gives £/ Cj qmc ~ 3.4i for 
the opening of a single-particle gap^ 7 -, while in the fRG 
we find U c « 2.6t as well (see Fig. [5]). 

The experimental study of Ref. [l9| does not find a gap 
for ABA(*) trilayer graphene. This is compatible with 
our findings, given the interactions in real material are 
weaker than this critical value. This is to be expected 
for consistency, as the critical interaction strength of the 
single layer and the ABA*-trilayer are roughly the same, 
and the single-layer remains semi-metallic, too. Even 
if they are slightly above the threshold, the expected 
transport gaps would be small and very hard to mea- 
sure. Therefore, while a precise quantitative picture can- 
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FIG. 6: Critical scale of ABA-trilayer including remote hop- 
pings (denoted ABA*, solid black line), with the ratios of the 
hopping parameters corresponding to the realistic value a 26 ' 27 , 
and with pure onsite repulsion U as interaction. For compar- 
ison, we also show the critical scales of ABA-trilayer without 
remote hoppings (dashed orange), AB bilayer (dot dashed 
blue), both with the choice t± = O.lt, and the single-layer 
system (dotted purple). The dashed horizontal line visualizes 
the energy scale A* where the bands become non-monotonous 
due to remote hoppings. 



not be obtained with our approximate fRG method and 
here for onsite interactions only, on a qualitative level 
we reach consistent conclusions. It would be interesting 
to resolve better the critical region close to C/ c .tl for the 
ABA* stacked trilayer with remote hoppings and analyze 
the onset of instabilities and their nature. However, this 
would require a different implementation of our patching 
scheme which we leave out for future work. In the re- 
mainder of this work, we will therefore concentrate on a 
more thorough study of the ABC trilayer model. 



B. ABC trilayer: Instabilities with non-local 
interactions 

From the results for pure onsite interactions and also 
inspired by recent experiments, we conclude that while 
ABC trilayers strongly support the formation of many- 
body states, the ABA* trilayer including remote hop- 
pings most probably does not. We take this as a motiva- 
tion to study further the instabilities of ABC trilayers for 
wider range of non-local interaction parameters, which 
brings us closer to the real materials. Ab initio values 
for the strength of the density-density interactions up 
to the third nearest neighbor were listed for single-layer 
graphene and graphite in Ref. [U Most likely, one can 
safely interpolate the parameters for the bi- and trilayer 
case from this data. 

Running the fRG for extended interactions, we find 
a number of different phases as previously described for 
the bilayer case in Ref. |32| . namely a charge density wave 
(CDW), a quantum spin Hall state (QSH) and charge 
density wave with non-zero momentum transfer (CDW3) 
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FIG. 7: Effective interaction vertex near the critical scale in 
the CDW regime for U = 0, Vi = 0.5t and V 2 = in units of 
t. Left Panel: Orbital combinations with 01 = 02 = 03=04. 
The numbers on the axis specify the number of the patch as 
shown in Fig. [3] On the horizontal axis the wavevector fci 
can be read off and on the vertical axis we enumerate ki- k% 
is fixed on the first patch, k$ then follows from momentum 
conservation. Left Panel: Effective vertex function for the 
orbital combination, where 01 = 03,02 = 04 ^ 01 and if 
01 G {01, aa, 03} then 02 G {61,62,63}- If in the right plot both 
01,02 G {(11,0,2,0,3} or G {61,62,63} the sign of the vertices 
changes. 

alongside the AF-SDW. For the investigation of the phase 
diagram we take into account non-local interaction con- 
tributions, namely the nearest-neighbor in-plane interac- 
tion Vi and the second-nearest neighbour in-plane inter- 
action V2 ■ More remote interaction contributions are ne- 
glected. In the bilayer case^ we checked explicitly that 
a third-nearest neighbor repulsion V3 does not change 
the picture. Also, we do not consider interlayer interac- 
tions. The ab-initio calculations in Ref. Hi] showed that 
the nearest interlayer interactions in graphite (as well as 
in layered graphene) are of the order of the V3-terrr>££. 

The two leading non-local interaction terms V\ and V2 
trigger the appearance of qualitatively different instabil- 
ities. For dominating V\ we find an instability toward a 
charge density wave with a momentum signature of the 
effective interaction as shown in Fig. [7] This momentum 
structure can be written down as an effective interaction 
Hamiltonian of the form, 

#cdw = -JfJ2 Voo>eo£ >N N°' (23) 
0,0' 

with V QO i > and N° = YV C* -c r. This sign struc- 

00 ^k,s 0:St k o,s,k o 

ture supports an enhanced occupancy of the Ai sublat- 
tices as compared to the Bi sublattices or vice versa. 
Furthermore, a mean-field decoupling of this effective in- 
teraction gives a gap in the single-particle spectrum. 

A dominating V2 yields an instability whose dominant 
interaction terms can be cast into an effective Hamilto- 
nian of the following type, 

Hqsh = ~ Yl V °°' e ° e °' &f ■ ®f ( 24 ) 

0,0' 

with V 00 , > and S° f = § Y.Z,,,* fk°™' c l >s ,k C o, s >,k 
including a /-wave form factor _/> = sin(v / 3afc a: ) — 
2sin(^%^)cos(^). This effective Hamiltonian can 



FIG. 8: Effective interaction vertex near the critical scale in 
the QSH regime for U = 0, Vi = and V2 = 1.5t in units of 
t. Left Panel: Orbital combinations with 01=02 = 03 = 04. 
The numbers on the axis specify the number of the patch as 
shown in Fig. [jj] On the horizontal axis the wavevector fei 
can be read off and on the vertical axis we enumerate k^. k% 
is fixed on the first patch, k$ then follows from momentum 
conservation. Here, we see that the sharp vertical structure 
[k\ = kz) comes with the double magnitude as the horizontal 
structure (k% = ^3). Middle Panel: Effective vertex function 
for the orbital combination, where 01 = 03,02 = 04 7^ 01 
and if 01 G {01,02,03} then 02 G {61,62,63}. Right panel: 
Effective vertex function for the orbital combination, where 
01 = 04,02 = 03/01 and if 01 G {01,02,03} then 02 G 
{61,62,63}. If in the second and third plot both 01,02 G 
{01,02, 03} or G {61, 62, 63} the sign of the vertices changes . 

be decoupled in a purely imaginary 'Kane-Mele' order 
parameter. This type of instability represents a many- 
body path to the quantum spin Hall state, where the 
mass term due to spin-orbit interaction in the original 
Kane-Mele proposal^ is now provided by an interaction- 
induced mean-field. In the single-layer— and bilayer— 
honeycomb models, this instability was found in the same 
corner of interaction parameter space. For an odd num- 
ber of layers, the mean-field Kane-Mele-ordered state will 
give rise to an odd number of helical edge modes, and 
thus represent a two-dimensional interaction-driven topo- 
logical insulator with protected edge modes. 

Finally, in the niche of the parameter space for smaller 
U, we also recover the CDW3 phase, that we already 
found in the bilayer system^, 

H C BW 3 ="^E V °°' e ° e °' ( N °Q N -Q + N °-Q N i) ( 25 ) 

O.O' 

with N% = E~ k ,s c l s r k+ Q C o,s,k- Scc Fig. El for the char- 
acteristic momentum structure of the effective interac- 
tion. The order parameter due to the symmetry break- 
ing (Nt) 7^ is in complete analogy to the one in the 

honeycomb bilayer—, except for the adapted definition 
of the e . Within one layer, this order forms three in- 
equivalent sites with different charge densities. The sign 
structure of the order parameter on different layers and 
sublattices is determined by the e -factors, so as to lower 
the energy contribution from ([23]). This leaves the total 
phase of the order parameter undetermined. Depend- 
ing on this phase, the quadratic mean-field Hamiltonian 
for a single layer either exhibits a gapless spectrum with 
Dirac points shifted away from the K, K' points or a 
fully gapped state. Which case represents the mean-field 
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FIG. 9: Effective interaction vertex near the critical scale in 
the CDW 3 regime for U = 0, Vi = and V 2 = 0.5i in units of 
t. Left Panel: Orbital combinations with 01=02 = 03=04. 
The numbers on the axis specify the number of the patch as 
shown in Fig. [jj] On the horizontal axis the wavevector fei 
can be read off and on the vertical axis we enumerate k^. 
ks is fixed on the first patch, k$ then follows from momen- 
tum conservation. Here, we see that the order parameter 
has non vanishing momentum transfer as the sharp feature 
is not at ki = kg. Right Panel: Effective vertex function for 
the orbital combination, where 01 = 03, 02 = 04 7^ 01 and if 
01 £ {0,1,0,2,0.3} then 02 £ {61,62,63}. If in the right plot 
both 01,02 £ {a\, 02, 03} or £ {61,62,63} the sign of the ver- 
tices changes. 



ground state as function of the interaction parameters 
has yet to be determined variationally. 



C. ABC trilayer phase diagram from fRG 

For a systematic investigation of the ABC trilayer 
phase diagram, we scan a range of interaction param- 
eters U, V\ and V2, whose ab-initio values are listed in 
Ref. As we expect the fRG to overestimate the crit- 
ical scales, we take these ab-initio parameters as upper 
bounds for the range our phase diagrams. In Fig. [101 wc 
show the fRG phase diagram obtained by identifying the 
leading tendencies in the effective interactions near the 
instability with an underlying contour plot of the critical 
scale A c in units of t. 

We also mark the ab-initio values in the lower plots of 
Fig. [TQl obtained by taking the values from Ref. |3l| and 
scaling {U, Vi, V2} — > &{U, V\, V2) so as to hit the values 
U = 2t and U = 3t. In both cases, these choices place the 
system near the phase boundary between QSH and AF- 
SDW instability. For the more long-ranged single-layer 
graphene interactions, one finds a QSH state, while for 
the slightly shorter ranged graphite parameters, one gets 
a AF-SDW state. Hence, which order occurs might be a 
delicate issue that is decided by details. In our approxi- 
mation, the critical scales interpolate smoothly across the 
phase borders, indicating a weaker competition between 
the different tendencies. Note that due to this borderline 
situation there is no true necessity for different layered 
graphene systems, e.g. with different environments, to 
exhibit the same ground state, and even the energy scales 
or gaps could come out similarly despite different states 
might be selected. Hence, regarding the leading insta- 
bility, the situation in the ABC trilayer is very similar 
to the one found previously in the Bcrnal-stackcd bilaycr 
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FIG. 10: Tentative fRG phase diagram of the ABC trilayer 
with t± = O.lt. The black lines are guides to the eye and 
seperate the different regimes. The contourplot encodes the 
critical scale at which the vertices diverge. The rescaled 
ab-initio interaction parameters for graphene (circles) and 
graphite (squares) of Ref. [3]] are shown in the lower plots. 



system. 

As mentioned before, we also find a rather exotic den- 
sity wave phase CDW3 with a tripling of the unit cell 
within the layers. This state is subject to further in- 
vestigation. However it occurs only at quite unrealistic 
corners of the parameter space with dominant non-local 
terms. Hence, we do not discuss it further here. 



D. Energy scales of the ordering phenomena 

In our previous study of the bilayer system^, we en- 
countered a problem of the current model studies that 
in principal also affects the present work. However, we 
will now also offer an explanation of what happens and 
how one should read the data in order to get reasonable 
agreement with and a more quantitative picture of the 
experiments. 

As shown in this work and also in the previous pa- 
per on the bilayer system^, the simple models employed 
by virtually all many-body approaches to interaction ef- 
fects in few-layer graphene can produce very large criti- 
cal scales. This can be already seen in Fig. [TQl where we 
indicate the values of the critical scales in units of the 
hopping t w 3 eV. These scales are huge for most choices 
of the interaction parameters. This is not surprising. In 
our and other theoretical approaches, the large scales are 
simply caused by the perfect particle-hole nesting of the 
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FIG. 11: FRG critical scale for ABC trilayer for the pure 
Hubbard model (blue solid) and rescaled ab initio parame- 
ters for graphite (red solid) and single-layer graphene (orange 
solid) as described in the text with t±_ = O.lt. For compari- 
son we also show the corresponding critical scale for the AB 
bilayer structure (dashed) and the singlelayer structure (dot- 
ted). For the Hubbard (blue curves) and the rescaled graphite 
parameters (red curves) we find the system to be in the AF- 
SDW phase, for the rescaled graphene parameters the system 
is in the QSH state (orange curves). The dashed horizontal 
line marks 10 meV, which is the order of magnitude, where 
the topology of the band structure in ABC trilayer graphene 
changes and remote hoppings become important. 



band structure. Furthermore, from the comparison with 
QMC calculations in the case of pure onsite interactions^ 
we know that the overestimate of the fRG is certainly 
not severe and, expressed conservatively, is less than an 
order of magnitude in the regime where also the QMC 
finds robust ordering. As the critical scale is an esti- 
mate for the energy scale of the spectral restructuring or 
gap opening in the ordered phase, a high critical scale 
would correspond to large energy gaps. If we took the 
ab-initio parameters of Ref. [3l| literally, the theoretical 
gap estimates would exceed the experimentally observed 
gap scales ~ 1 — 10 meV by orders of magnitude. 

Let us now analyze the systematics of these critical 
scales a bit more deeply. In Fig. [TTJ we show the critical 
scales obtained from fRG for the single-, bi- and trilayer 
graphene for three cases, namely for onsite interaction U 
only, and for nonlocal interactions with the cRPA param- 
eters V\ and V2 for graphite and graphene with repulsions 
up to the second nearest neighbor, for realistic interlayer 
hopping t± = O.lt. The curves show the dependence 
on the overall magnitude of the onsite interaction, where 
the ratio between the local U and the non-local inter- 
action parameters is kept fixed. Obviously there are two 
regimes: a high-scale regime with large critical scales that 
do not depend too strongly on the interaction strength 
and a low-scale regime with a very strong dependence on 
the interaction. For the single layer the second regime is 
very narrow and basically only contain the minimal criti- 
cal interaction strength U c (Vi, V2) below which the semi- 
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FIG. 12: FRG critical scale for ABC trilayer (solid lines), 
ABA trilayer without remote hoppings (dashed lines) and AB 
bilayer (dotted lines) with fixed onsite interaction and variable 
t±. The red lines show the results for U = 3t > U c ,sh and no 
dependence of the critical scale on the interlayer hopping t± ■ 
For the case U — It < U c ,sl (black lines) we observe a strong 
dependence of the critical scale on t±. 



metal is stable. Also for bi- and trilayer, the high-scale 
regime starts above the single-layer U c (Vi, V2). Next let 
us consider the dependence of the critical scale A c on the 
perpendicular hopping tj_ in the ABC trilayer. Whereas 
the nature of the ground state qualitatively remains the 
same for all choices of t± ^ the absolute value of this 
parameter has an impact on the size of the critical scale 
and shows different behaviors on the two different sides 
of the critical onsite interaction of single-layer graphene 
U c .sL, sec Fig.Q/Jl i.e. whether we are in the high-scale or 
in the low-scale regime. For large interactions U > U c .sl 
the size of t± is of no importance for the critical scale A c . 
This changes for smaller U. Here, the smaller t± is, the 
stronger is the A c -variation with the interaction strength. 
The comparison with QMC in Ref. [33| was done at larger 
t± = t and U > 2.8t, where the scales do not vary that 
strongly. For band structure parameters with ij_ < O.lt 
and small interactions U < U Ci sl, we also observe flows 
to strong coupling for the ABC and the ABA trilayer 
system, however, the critical scales turn out to be very 
small, an effect that reflects the behavior of the single- 
layer system, where no instabilities occur for U < U Cj sl- 
While this does not constitute a difficulty for the fRG 
method per se it makes the numerical evaluation very 
tedious and renders the comparison with QMC impossi- 
ble. 

We now argue that in order to account for the approx- 
imations made in the fRG scheme, mainly the neglect of 
self-energy corrections, and in order to obtain a realistic 
picture, we have to reduce the cRPA parameters with 
U ~ 3t by hand. The reduction factor is chosen so as 
to obtain the experimentally verified semi-metallic solu- 
tion for the single-layer case, i.e. roughly to U ~ 1.5t 
(if we take the single-layer cRPA values) with an anal- 
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ogous rescaling of the non-local couplings by the same 
factor. This shift on the interaction axis now takes us 
from the high scale regime into the regime of strongly 
varying scales, cf. Fig. [TT] Now the critical scales for 
bi- and trilayer end up in the range < O.Olt ~ 30meV 
which is much closer to the experimental values for gap 
sizes and already in the correct order of magnitude. 

In the experiments on bi- and trilayer graphene of these 
types, the energy gap in the trilayer— came out a factor 
2-3 higher than in the bilayer. In the fRG the ratio be- 
tween the critical scales for these two systems depends on 
the parameter values for the interactions. However for in- 
tercations where single-layer graphene does not undergo 
a phase transition, the larger density of states near the 
Fermi level of the trilayer case leads to a larger critical 
scale in theory as well, see Fig. [TTJ Hence the different 
energy gaps in bi- and trilayers are qualitatively captured 
by this theory. Remarkably, with this choice of interac- 
tion parameters, the instabilities in ABC trilayer occur 
on an energy scale where remote hoppings start to be 
important (~ lOmcV, see Fig. [Til. 

V. DISCUSSION 

We have performed extensive fRG calculations on hon- 
eycomb trilayer systems with different stacking orders, 
as model systems for trilayer graphene. In doing this 
we have used as far as possible the available input pa- 
rameters from ab-initio calculations. Moreover, we have 
taken into account the full 6-band band structure ob- 
tained within the model with one effective p 2 -type orbital 
per carbon site. The fRG approach we have used is cer- 
tainly not exact, but it goes far beyond mean-field stud- 
ies, random-phase approximation approaches and other 
perturbative calculations performed for these systems. 
Furthermore, in the bilayer case with onsite interactions 
only, the qualitative and quantitative information of the 
fRG compares quite favorably with the results of QMC 
calculations^. For non-local interactions, QMC encoun- 
ters severe sign problems, and there is no possibility for 
benchmarks, but there is also no clear reason why the 
fRG method should get worse then. Hence, for the given 
theoretical model the fRG results should be rather ro- 
bust both qualitatively and also quantitatively (with the 
adaptations discussed is Sec. HVDI and mentioned below) 
at least as an order-of-magnitudc estimate. 

First let us mention the main results of the present 
trilayer study before we get to the connections with 
graphene systems and experiments. Comparing the dif- 
ferent trilayer stackings with the Bernal-stacked bilayer 
and the monolayer results, we could identify the ABC 
trilayer as the system that is most prone toward insta- 
bilities, with larger energy scales for ordering than the 
AB-bilayer. The ABA trilayer without the additional in- 
terlayer 'remote' hoppings is comparable in its critical 
scales to the bilayer, but we showed the remote hop- 
ping terms with the suggested realistic parameter val- 



ues (called ABA* structure here) are likely to remove 
the instability at least for smaller, possibly realistic, in- 
teraction strengths. Interestingly, for the Hubbard on- 
site interaction case, the minimal £/-value for obtaining 
a gapped ground state in the ABA*-structure is close to 
the one for the single layer. Taking the current exper- 
imental knowledge for bi- and trilayer graphene, these 
theoretical findings regarding the systematic differences 
are fully consistent with the observations. 

Due to the uncertainties about the parameter values 
for the theoretical model and the approximations made 
in our approach, we cannot expect a fully quantitative 
description. However, the phenomenological input of 
requiring the single-layer to remain semi-metallic puts 
bounds on the bare interaction parameters that should 
be used in our model. The idea used in this paper is 
to scale down the ab-initio parameters for the interac- 
tions by an appropriate factor in order to compensate 
for the approximations, so as to render the fRG flow for 
the single layer free of divergences. With this remedy 
for the inexactness of our approach, the ABA* trilayer 
with remote hoppings remains semi-metallic and the en- 
ergy scales of the unstable systems, namely the bilayer 
and the ABC trilayer come out in a quite realistic region 
below 30 meV. 

The type of the leading instability in these trilayer 
systems can also be read off from the fRG effective in- 
teractions near the instability. Here, for the interaction 
parameters determined by ab-initio methods^, our cal- 
culations show a strong competition between SDW anti- 
ferromagnetic ordering and the interaction-driven 'Kane- 
Mele' quantum spin Hall state. Which tendency wins 
depends on the detailed spatial profile of the interac- 
tions. A more longer-ranged behavior favors the QSH 
instability, and for pure onsite repulsions, the AF-SDW 
state is the clear winner. Both states would open up a 
bulk gap. The SDW states should have an interesting 
modulation of the ordered moments depending on the 
number of nearest neighbors, with smaller moments for 
higher coordination number—. The QSH state should be 
a true two-dimensional topological insulator for the tri- 
layer case, as time-reversal invariant edge defects will not 
suffice to gap the three pairs of helical edge states com- 
pletely, in contrast with the bilayer case (for a discussion 
of the edge-state robustness, see the review in Ref. [50h . 
Hence, at least one pair of counter moving, spin-resolved 
helical edge states should survive time-revcrsal-invariant 
edge disorder and could hence serve as smoking gun for 
such a correlated ground state. This perspective is rather 
exciting, as this would be the first realization of the con- 
cept of an interaction-driven topological insulator ('topo- 
logical Mott insulator')^. 
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